Cite this Article
Huang Li. Kernel polynomial representation for imaginary-time Green’s functions in continuous-time quantum Monte Carlo impurity solver. Chinese Physics B, 2016, 25(11): 117101  
Permissions
Kernel polynomial representation for imaginary-time Green’s functions in continuous-time quantum Monte Carlo impurity solver
Huang Li†,
Science and Technology on Surface Physics and Chemistry Laboratory, China Academy of Engineering Physics, Jiangyou 621908, China

 

† Corresponding author. E-mail: lihuang.dmft@gmail.com

Project supported by the National Natural Science Foundation of China (Grant No. 11504340).

Abstract
Abstract

Inspired by the recently proposed Legendre orthogonal polynomial representation for imaginary-time Green’s functions G(τ), we develop an alternate and superior representation for G(τ) and implement it in the hybridization expansion continuous-time quantum Monte Carlo impurity solver. This representation is based on the kernel polynomial method, which introduces some integral kernel functions to filter the numerical fluctuations caused by the explicit truncations of polynomial expansion series and can improve the computational precision significantly. As an illustration of the new representation, we re-examine the imaginary-time Green’s functions of the single-band Hubbard model in the framework of dynamical mean-field theory. The calculated results suggest that with carefully chosen integral kernel functions, whether the system is metallic or insulating, the Gibbs oscillations found in the previous Legendre orthogonal polynomial representation have been vastly suppressed and remarkable corrections to the measured Green’s functions have been obtained.

PACS: 71.10.Fd;71.27.+a
Keyword:kernel polynomial representation;imaginary-time Green’s function;continuous-time quantum Monte Carlo impurity solver;dynamical mean-field theory
1. Introduction

The rapid development of efficient numerical and analytical methods for solving quantum impurity models has been driven in recent years by the unprecedented success of dynamical mean-field theory (DMFT)[1–3] and its non-local extensions.[4–6] In the framework of DMFT, the self-energy function is assumed to be local; i.e, its momentum dependence is completely neglected. Thus, the solution of the original lattice model may be obtained from the solution of an appropriately defined quantum impurity model plus a self-consistency condition. In order to solve the quantum impurity models, numerous quantum impurity solvers have been developed in the last few decades.[1,2] In particular, the recently developed continuous-time quantum Monte Carlo impurity solvers[7–12] are being more generally preferred due to their accuracy, efficiency, and ability to treat extreme low temperature and arbitrary interaction terms in the quantum impurity models.

Among various continuous-time quantum Monte Carlo algorithms, the hybridization expansion version (abbreviated CT-HYB)[8–10] is up to now the most powerful and reliable impurity solver and it is widely used. However, for the CT-HYB quantum impurity solver, several severe technical limits still remain. One well-known problem is the high frequency noise commonly observed in the Matsubara Green’s function G(iωn) and the self-energy function Σ(iωn).[13,14] Similar problems also arise in the calculations of imaginary-time Green’s functions G(τ) and two-particle vertex functions Γ(ω, ω′, ν),[15] which are less emphasized in the literature, according to our knowledge. In order to cure these problems, an intuitive and straightforward idea is to run more statistics in the Monte Carlo simulations to suppress the data fluctuations as far as possible. Actually, although this strategy mitigates the problems, it will not avert them and will deteriorate the efficiency of the CT-HYB quantum impurity solver rapidly.

Recently, Lewin Boehnke et al.[16] suggested to measure the single-particle and the two-particle Green’s functions in the basis of Legendre orthogonal polynomials. The Legendre orthogonal polynomial method provides a much more compact representation of Green’s functions than standard Matsubara frequencies and, therefore, significantly reduces the memory requirement for these quantities. Moreover, it can be used as an efficient noise filter for various physical quantities within the CT-HYB quantum impurity solver since the statistical noise is mostly carried by the high-order Legendre coefficients which should be truncated, while the physical properties are largely determined by the low-order coefficients which should be retained. By using the Legendre orthogonal polynomial representation, the accuracy of the CT-HYB quantum impurity solver is greatly improved. This representation has been broadly used in the computations of the single-particle Green’s function and lattice susceptibilities in the context of realistic DMFT calculations in combination with the density functional theory.[17,18]

It seems that the Legendre orthogonal polynomial representation is superior to the old one. However, according to careful examination, we find that the so-called Gibbs oscillations appear quite easily in the calculated Green’s functions and other physical quantities, which could be due to the rough truncation of the Legendre basis. The character of the Gibbs oscillations is that the calculated physical quantities are smooth but oscillating periodically with the scattered direct measurements. The situation is even worse when the system is in the insulating state where the Gibbs oscillations will likely cause the reconstructed single-particle Green’s function to break the causality. Note that a common procedure to damp these oscillations relies on an appropriate modification of the expansion coefficients by some integral kernel functions, which is the well-known kernel polynomial representation.[19] Motivated by this idea, in this study we adopt the kernel polynomial method (dubbed KPM) to improve the measurement of the single-particle and two-particle quantities for the CT-HYB quantum impurity solver. We find that significant improvements are gained, and the kernel polynomial representation is better than the Legendre orthogonal polynomial representation, no matter whether the system is metallic or insulating.

The rest of this paper is organized as follows. In Subsection 2.1, a brief introduction to the orthogonal polynomial representation is provided. The original implementation is based on the Legendre polynomials only. Here we generalize it to the Chebyshev polynomials. In Subsection 2.2, the kernel polynomial representation is presented. In Section 3, we benchmark the KPM by re-examining the imaginary-time Green’s functions G(τ) for the single-band Hubbard model with half-filling. The characteristics of different integral kernel functions which are used to adjust the expansion coefficients are discussed in detail. Section 4 gives a conclusion and outlook.

2. Method
2.1. Chebyshev orthogonal polynomial representation

In principle, the imaginary-time Green’s function G(τ) (τ ∈ [0, β]) can be expanded in terms of orthogonal polynomials defined on the interval [−1, 1], such as the Legendre polynomials Pn(x),

where β is the inverse temperature (≡ 1/T), x(τ) = 2τ/β − 1, and Gn denotes the expansion coefficients of G(τ) in the Legendre orthogonal polynomials basis[16]

Since the expansion coefficients Gn generally show a very fast decay with n, the expansion in the Legendre polynomials can be truncated at a maximum order nmax. In the CT-HYB quantum impurity solver, the formula for measuring imaginary-time Green’s function G(τ)[8–10] reads

where k is the order of diagrammatic perturbation expansion series, matrix element with Δ(τ) being the hybridization function, and and are the coordinates in the imaginary-time axis for creation and annihilation operators, respectively. By utilizing Eqs. (3) and (4), the Legendre expansion coefficients for G(τ) finally become

Lewin Boehnke et al.[16] have chosen the Legendre orthogonal polynomials as their preferred basis to expand the single-particle and the two-particle Green’s functions. But it should be emphasized that a priori different orthogonal polynomial basis is possible as well. Thus we would like to generalize the approach to use the Chebyshev orthogonal polynomials as an optional basis. It is well-known that there exist two kinds of Chebyshev polynomials.[20] Here we choose the second kind Chebyshev polynomials Un(x) as the basis because their orthogonal condition is simpler. Then, we reach the following expressions:

After a straightforward substitute, the Chebyshev expansion coefficients for G(τ) finally become

where

Virtually, the CT-HYB quantum impurity solver can directly accumulate the Legendre or Chebyshev expansion coefficients Gn instead of the original Green’s functions G(τ). Once the Monte Carlo sampling has been finished, G(τ) can be reconstructed analytically via the expansion coefficients (see Eqs. (1) and (7)). Since the expansion coefficients decay very quickly, only a small number of coefficients are needed. As a result, the orthogonal polynomial bases are much more compact and particularly useful for storing and manipulating the two-particle quantities, such as the two-particle vertex function, etc. Furthermore, the Monte Carlo noises are mainly concentrated in the high-order expansion coefficients, and their numerical values are usually very small. So a rough truncation method can be developed to filter out the noises and obtain more smooth and accurate results. The performances of Legendre and Chebyshev polynomials are similar. In the next section, we will prove this using some concrete examples.

2.2. Kernel polynomial representation

As we already know, the essential idea of the orthogonal polynomial method is to expand G(τ) and the other physical observables in infinite series of Chebyshev or Legendre polynomials, and then use the Monte Carlo algorithm to sample the expansion coefficients Gn directly. However, the infinite expansion series will remain finite actually from a numerical perspective (n ≤ nmax), and we are thus trapped by a classical problem of approximation theory. In this case, the problem is equivalent to figure out the best approximation to G(τ) given a finite number of Gn. Experience shows that a crude truncation of the infinite series probably leads to poor numerical precision and fluctuations, which is also known as Gibbs oscillations in the literature.[19] Later, we will see that as for the reconstructed Green’s function G(τ) in the insulating state, almost periodic Gibbs oscillations are clearly identified in a wide τ range.

The KPM is an efficient approach to systematically damp the Gibbs oscillations. Its spirit is to introduce some kind of integral kernel function fn, and then the expansion coefficients are modified from Gn to Gnfn.[19] The fn decreases monotonously with respect to n and satisfies 0.0 ≤ fn ≤ 1.0, so as to Gn fn → 0.0 when n → ∞. In order words, fn can be considered as a damping factor. Obviously, the simplest integral kernel function, which is usually attributed to Dirichlet, is obtained by setting fn ≡ 1. By using the Dirichlet kernel, the KPM is equivalent to the previous orthogonal polynomial method. In addition to the Dirichlet kernel, there also exist many other integral kernel functions, such as Jackson, Lorentz, Fejér, and Wang–Zunger kernels. Their formulas are collected in Table 1 and plotted in Fig. 1, respectively. Note that for all of the kernel functions, f0 must be equal to 1 and fn must approach 1 as n → ∞. The optimal integral kernel function partially depends on the considered problem. According to the literature,[19] the Jackson kernel may be the best for most applications, and the Lorentz kernel may be the best for Green’s function.

Fig. 1. Typical integral kernel functions fn used to improve the quality of the polynomial expansion series. The order of expansion series is fixed to N = 64. The Dirichlet, Jackson, and Fejér kernels take no parameters. For the Lorentz kernel, the λ parameter is fixed at 1.0. For the Wang–Zunger kernel, parameters α and β are 1.0 and 4.0, respectively.

Finally, we note that the integral kernel functions fn can be evaluated and stored in advance, so that the KPM has no effect on the computational efficiency of the CT-HYB quantum impurity solver. The implementation of KPM is very simple, only trivial modifications are needed for the orthogonal polynomial representation version of CT-HYB quantum impurity solver, i.e., replacing Gn by Gnfn. Since both the kernel polynomial and orthogonal polynomial representations are only alternate bases for the single-particle and two-particle quantities, they can be implemented in the segment picture[8] and general matrix[9] formulations of CT-HYB quantum impurity solvers to improve the numerical accuracy and efficiency. We have implemented the kernel polynomial and orthogonal polynomial representations (using the Legendre and Chebyshev orthogonal polynomials) in a generic quantum impurity solver package iQIST.[21]

Table 1.

Collection of various integral kernel functions fn that can be used to improve the quality of an order-N Legendre or Chebyshev series.

.
3. Benchmark

In this section we will try to benchmark the kernel polynomial representation and compare the calculated results with those obtained by the orthogonal polynomial representation[16] and conventionally direct measurements.[8] For the sake of simplicity, a single-band Hubbard model with half-filling is chosen as a toy model to examine our implementations. The Hamiltonian of the single-band Hubbard model reads

where is the local Hamiltonian on each lattice site i

Here σ is the spin index, t is the hopping parameter, μ is the chemical potential, and U is the intra-orbital Coulomb interaction. In the present study, we just consider two typical scenarios: (i) U = 4.0 and β = 10.0 for the metallic case, and (ii) U = 6.0 and β = 50.0 for the insulating case. The bandwidth is 2.0 (= 4t), and a semicircular density of states is chosen. The chemical potential μ is fixed at U/2 to meet the half-filling condition. This model is studied within the framework of single site DMFT,[1,2] and the segment picture version of CT-HYB[8] as implemented in the iQIST software package[21] is used as the quantum impurity solver. In each DMFT iteration, typically 4 × 108 Monte Carlo sweeps are performed to reach a sufficient numerical precision. We also carry out additional ED calculations to generate accurate imaginary-time Green’s function as a reference.

3.1. Metallic state

Let us first concentrate our attention on the metallic state. In Fig. 2, the “bare” expansion coefficients β|Gn| for the Legendre and Chebyshev orthogonal polynomials are shown. The characteristics of the expansion coefficients for the two kinds of orthogonal polynomials are very similar, which means that the performance (convergence behavior) of the Chebyshev polynomial representation is comparable to that of the Legendre polynomial representation. Next we analyze their common features. As is already pointed out by Lewin Boehnke et al.,[16] due to the constraint of the particle–hole symmetry, the expansion coefficients for odd order n should be zero. Indeed, the coefficients in our data for odd n all take very small values, compatible with a vanishing value within their error bars. The even n coefficients instead show a very fast decay. As is shown in this figure, nmax = 15–25 is enough for both the Legendre and Chebyshev polynomial representations to obtain converged and accurate results. Thus, in this case, nmax is fixed to be 24.

Fig. 2. Calculated Legendre and Chebyshev coefficients β|Gn| for imaginary-time Green’s functions of the single-band half-filled Hubbard model on the Bethe lattice within DMFT. The Coulomb interaction strength U is 4.0 and the inverse temperature β is 10. The coefficients in the pink region contribute very little to the resulting Green’s function.
Fig. 3. The imaginary-time Green’s function G(τ) in τ ∈ [0.4β, 0.6β] interval of the single-band half-filled Hubbard model. The Coulomb interaction strength U is 4.0 and β = 10. (a) G(τ) calculated by Legendre polynomials with or without integral kernel functions. (b) G(τ) calculated by Chebyshev polynomials with or without integral kernel functions. The maximum order for the polynomial expansion series, nmax, is 24. The ED data are used for comparison. All data are magnified by a factor of 100.

Figure 3 shows the calculated imaginary-time Green’s function G(τ) by using KPM with different orthogonal polynomials and integral kernel functions. It is apparent that the directly measured G(τ) is full of noise and fluctuation, which are bad for the later analytical continuation procedure.[22] Once the orthogonal polynomial method is used (i.e., the Dirichlet kernel fn = 1 is adopted), G(τ) turns smooth but obvious undulations still exist. If the Lorentz and Fejér kernels are applied respectively, the Green’s functions are smooth and without obvious undulations but deviate systematically from the directly measured data. If the Wang–Zunger kernel is used, the calculated G(τ) still exhibits a sizable undulation. Since the Wang–Zunger kernel is close to the Dirichlet kernel when n ≤ 20 (see Fig. 1), it is not surprising that the Green’s functions obtained by using these two kernels show the same oscillating behaviors. In a general view, the Jackson kernel function is the optimal choice. The resulting G(τ) evaluated by the Jackson kernel function is smooth and nicely interpolates the directly measured data. The Jackson kernel function has another important advantage in that it is parameter-free. On the contrary, the Lorentz and Wang–Zunger kernel functions need addition parameters (see Table 1). As is expected, the type of orthogonal polynomials has little impact on the interpolated G(τ). It seems that the Chebyshev polynomials do a bit better than the Legendre polynomials.

3.2. Insulating state

Now let us turn to the insulating state. When U = 6.0 and β = 50, a definitely insulating solution is obtained within DMFT. The “bare” Legendre and Chebyshev coefficients β|Gn| of G(τ) are shown in Fig. 4. Similar to the metallic state, Gn takes very small value for odd n and can be ignored safely, and for even n, Gn converges to zero very quickly. For the Chebyshev and Legendre polynomials, nmax = 35–45 and nmax = 40–50, respectively. Thus the Chebyshev polynomials are more compact and efficient than the Legendre polynomials. In this case, nmax is fixed to be 96 uniformly, which is enough to obtain smooth and accurate results.

Fig. 4. Calculated Legendre and Chebyshev coefficients β|Gn| for imaginary-time Green’s functions of the single-band half-filled Hubbard model on the Bethe lattice within DMFT. The Coulomb interaction strength U is 6.0 and the inverse temperature β is 50. The coefficients in the pink region contribute very little to the resulting Green’s function.
Fig. 5. The imaginary-time Green’s function G(τ) of the single-band half-filled Hubbard model. The Coulomb interaction strength U is 6.0 and β = 50. (a) G(τ) calculated by Legendre polynomials with or without integral kernel functions. (b) G(τ) calculated by Chebyshev polynomials with or without integral kernel functions. In panels (c) and (d), the fine structures of G(τ) in τ ∈ [0.2β, 0.8β] interval are shown. They are obtained by Legendre polynomials and Chebyshev polynomials, respectively. The maximum order for the polynomial expansion series, nmax, is 96. The ED data are used for comparison. All data in panels (c) and (d) are magnified by a factor of 10000 for a clear view.

The calculated imaginary-time Green’s function G(τ) by using KPM with different orthogonal polynomials and integral kernel functions are illustrated in Fig. 5. A cursory look could lead us to believe that the reconstructed Green’s functions by the KPM or orthogonal polynomial method agree rather well with the directly measured data. However, this intuition is not correct. Next let us zoom in τ ∈ [0.2β, 0.8β] interval and check carefully the magnified G(τ), which are just depicted in Figs. 5(c) and 5(d). Clearly, G(τ) is approximately flat and very close to zero in this region. The default data obtained by the direct measurement exhibit periodical undulations. The reconstructed G(τ) by the orthogonal polynomial method (using the Dirichlet kernel) displays stronger periodical oscillations and violates the causality at the same time, which means that the results will be even deteriorated by using the orthogonal polynomial representation. The results obtained by the Wang–Zunger kernel fit the default data quite well, but obvious improvement is not observed. For the Lorentz and Fejér kernels, the calculated results under-fit the Green’s function and deviate from the reference data systematically. Again, the Jackson kernel function is the optimal choice. The calculated results are quite smooth, perfectly interpolate the reference data, and obey the causality all the time.

4. Conclusion

It is suggested that the orthogonal polynomial method based on the Legendre orthogonal polynomials can be used to improve the accuracy and computational efficiency of the CT-HYB quantum impurity solver,[16] but it still suffers from the Gibbs oscillations since the number of polynomial expansion coefficients is limited. In this paper, we develop a better representation to cure this problem.

At first, we generalize the original orthogonal polynomial method to using the second kind Chebyshev orthogonal polynomials. We find that the Chebyshev polynomials require fewer expansion coefficients to reach the converged results, and thus are better than the Legendre polynomials in computation efficiency. Second, we believe that the KPM with carefully chosen integral kernel functions could be used to suppress the Gibbs oscillations observed in the single-particle Green’s functions obtained by the orthogonal polynomial method and improve the numerical accuracy further. According to the benchmark results for the single-band half-filled Hubbard model, it is demonstrated that the Jackson kernel function is the optimal choice for imaginary-time Green’s function G(τ). Finally, though the KPM presented in this paper is mainly developed for the CT-HYB quantum impurity solver, it can be easily generalized to other continuous-time quantum Monte Carlo impurity solvers.[7] Besides the imaginary-time Green’s function, this method can be applied to the other physical observables, such as Matsubara Green’s function G(iωn) and self-energy function Σ(iωn), the generalization is straightforward.

Reference
1Georges AKotliar GKrauth WRozenberg M J 1996 Rev. Mod. Phys. 68 13
2Kotliar GSavrasov S YHaule KOudovenko V SParcollet OMarianetti C A 2006 Rev. Mod. Phys. 78 865
3Held K 2007 Adv. Phys. 56 829
4Maier TJarrell MPruschke THettler M H 2005 Rev. Mod. Phys. 77 1027
5Toschi AKatanin A AHeld K 2007 Phys. Rev. 75 045118
6Rubtsov A NKatsnelson M ILichtenstein A I 2008 Phys. Rev. 77 033101
7Emanuel GMillis A JLichtenstein A IRubtsov A NTroyer MWerner P 2011 Rev. Mod. Phys. 83 349
8Werner PComanac Ade’Medici LTroyer MMillis A J 2006 Phys. Rev. Lett. 97 076405
9Werner PMillis A J 2006 Phys. Rev. 74 155107
10Haule K 2007 Phys. Rev. 75 155113
11Gull EWerner PParcollet OTroyer M 2008 Europhys. Lett. 82 57003
12Rubtsov A NSavkin V VLichtenstein A I 2005 Phys. Rev. 72 035122
13Gull EWerner PMillis A JTroyer M 2007 Phys. Rev. 76 235123
14Blumer N 2007 Phys. Rev. 76 205120
15Kune Å J 2011 Phys. Rev. 83 085102
16Boehnke LHafermann HFerrero MLechermann FParcollet O 2011 Phys. Rev. 84 075145
17Deng X YFerrero MMravlje JAichhorn MGeorges A 2012 Phys. Rev. 85 125137
18Boehnke LLechermann F 2012 Phys. Rev. 85 115128
19Weiße AWellein GAlvermann AFehske H 2006 Rev. Mod. Phys. 78 275
20Gradshteyn I SRyzhik I MJeffrey AZwillinger D2000Table of Integrals, Series, and Products6th EdnAmsterdam Academic
21Huang LWang Y LMeng Z YDu LWerner PDai X 2015 Computer Physics Communications 195 140
22Jarrell MGubernatis J E 1996 Phys. Rep. 269 133